#### Conceptualizing and Measuring Support for Democracy: A New Approach
#### R code to run cross-country analyses using national sample results

library(psych)
library(lme4)
library(reshape2)
library(tidyverse)
library(pals)
library(countrycode)
library(colorspace)
library(XML)
library(RColorBrewer)

# options
par(mfrow=c(1,1), mar=c(3, 3, 0.5, 0.5), tcl=-0.2, las=0, mgp=c(1.5, 0.6, 0))

# country names
cnts = c("Argentina", "Brazil", "Britain", "Chile", "Germany", "Greece", "Hungary", "Israel", 
         "Mexico", "Netherlands", "Norway", "Peru", "Poland", "Portugal", "South Africa", 
         "Spain", "Taiwan", "Turkey", "USA")
n.cnts = length(cnts)



#### Plot 17-item CFA loadings

CFA_tables = list()

# Store data from each country's results in a list
for (i in 1:length(cnts)) {
  file = paste0("./", cnts[i], "/cfa_1fac_ord_std_table.html")
  # Attempt to read the HTML table and handle potential errors
  tryCatch({
    tab_temp = readHTMLTable(file, which = 1)
    tab_temp = tab_temp[3:20, c("V1", "V2")]
    # rename the column V2 with the country name
    colnames(tab_temp)[2] = cnts[i]
    CFA_tables[[i]] = tab_temp
  }, error = function(e) {
    cat("Error for", cnts[i], ": ", e$message, "\n")
  })
}

# merge, rearrange rows, and clean up
CFA_merg = Reduce(function(x, y) merge(x, y, by = "V1", all = TRUE), CFA_tables)
row_lab = c("FREXP1", "FREXP2", "FRASSC1", "FRASSC2", "FRASSC3", "UNISUFF1", "UNISUFF2", "DECELEC1", 
              "DECELEC2", "FRELECT1", "FRELECT2", "JUDCNSTR1", "JUDCNSTR2", "LEGCNSTR1", "LEGCNSTR2", 
              "EQLAW1", "EQLAW2")
CFA_merg$V1 = factor(CFA_merg$V1, levels = row_lab)
CFA_merg = arrange(CFA_merg, V1)
CFA_merg = filter(CFA_merg, V1 != "")
rownames(CFA_merg) = CFA_merg[,1]
CFA_merg[,1] = NULL


## Set up plot options

n.cnts = length(cnts)

# country code
cnts1 = cnts
cnts_code = countrycode(cnts1, origin="country.name", destination="iso2c")

# title-case country names
cnt_plot = str_to_title(cnts1)
cnt_plot[n.cnts] = "USA"

# more meaningful labels
it_labs = c("FreeExp1(r)", "FreeExp2", "FreeAssc1", "FreeAssc2(r)", "FreeAssc3", "UniSuff1", 
            "UniSuff2(r)", "ElecDecMk1", "ElecDecMk2(r)", "FFElect1(r)", "FFElect2", 
            "JudCnstr1(r)", "JudCnstr2", "LegCnstr1", "LegCnstr2(r)", "EqLaw1", "EqLaw2(r)")

# plot colours
plot_cols = unname(alphabet(n=n.cnts))
plot_cols[4] = "#740AFF"

# horiz lines separating components
comp_bounds = c(2.5, 4.5, 6.5, 8.5, 10.5, 12.5, 15.5)


## Plot using component labels

comp_labs = c("Freedom of expression", "Freedom of association", "Universal suffrage", 
              "Key decision-makers elected", "Free and fair elections", "Judicial constraints",
              "Legislative constraints", "Equality before the law")

cfa_dat_2 = matrix(NA, nrow=25, ncol=19)
lab_inds = c(1, 4, 8, 11, 14, 17, 20, 23)
dat_inds = c(2, 3, 5, 6, 7, 9, 10, 12, 13, 15, 16, 18, 19, 21, 22, 24, 25)
comp_pos = c(25, 22, 18, 15, 12, 9, 6, 3)
lab_pos = c(24, 23, 21, 20, 19, 17, 16, 14, 13, 11, 10, 8, 7, 5, 4, 2, 1)
for(i in 1:19) {
  cfa_dat_2[lab_inds, i] = NA
  cfa_dat_2[dat_inds, i] = CFA_merg[, i] 
}
it_labs_2 = c("Criticize govt.(r)", "Free media", "One party", "Protest(r)", "Ban orgs", "Voter competence", 
              "Tolerance(r)", "Technocratic rule", "Unelected authority(r)", "Respect results(r)", 
              "Bend rules", "Judicial review(r)", "Judicial compliance", "Legisl. compliance", 
              "Legis. review(r)", "Law enforced", "Access justice(r)")
comp_headers = lab_inds

pdf("fig_2.pdf", height=7, width=8)
par(mfrow=c(1,1), mar=c(2.5, 7.2, 0.5, 1), tcl=-0.2)
plot(x=cfa_dat_2[, 1], y=25:1, type="n", pch=16, col=rgb(0, 0, 0, 0.5), 
     xlim=c(-0.45, 1), xaxs="i", axes=FALSE, xlab="", ylab="")
axis(side=1, at=round(seq(-0.6, 0.8, by=0.2), 1), mgp=c(1, 0.2, 0), cex.axis=0.8)
axis(side=2, at=lab_pos, labels=it_labs_2, las=1, mgp=c(1, 0.4, 0), cex.axis=0.8)
mtext(side=1, line=1.2, text="Liberal democracy factor loading")
abline(h=lab_pos, lty=3, lwd=0.75, col=rgb(0.3, 0.3, 0.3, 0.75))
abline(v=seq(-0.6, 0.8, by=0.2), lty=3, lwd=0.75, col=rgb(0.3, 0.3, 0.3, 0.75))
abline(v=0, lty=1, lwd=0.75, col=rgb(0.3, 0.3, 0.3, 0.75))
abline(h=comp_pos, lty=1, lwd=25, col=rgb(0.7, 0.7, 0.7, 1))
text(x=0.275, y=comp_pos, adj=0.5, cex=0.9, labels=comp_labs)
for(i in 1:n.cnts){
  text(x=cfa_dat_2[, i], y=jitter(25:1, factor=1), labels=cnts_code[i], col=rgb(0, 0, 0, 1), cex=0.65)
}
box()
dev.off()



### 17-item CFA fit plot

## 1 factor fits

fit_list_1f = list()

# Store data from each country's results in a list
for (cnt in cnts) {
  file_path = paste0(cnt, "/cfa_1fac_std_ord_fit.csv")
  fit_dat_1f = read.csv(file_path, row.names = 1)
  fit_list_1f[[cnt]] = fit_dat_1f
}

# Combine the data from all files into a matrix
fit_mat_1f = do.call(cbind, fit_list_1f)
names(fit_mat_1f) = cnts
fit_mat_1f


## load 2 factor fits

fit_list_2f = list()

# Store data from each country's results in a list
for (cnt in cnts) {
  file_path = paste0(cnt, "/cfa_2fac_std_ord_fit.csv")
  fit_dat_2f = read.csv(file_path, row.names = 1)
  fit_list_2f[[cnt]] = fit_dat_2f
}

# Combine the data from all files into a matrix
fit_mat_2f = do.call(cbind, fit_list_2f)
names(fit_mat_2f) = cnts
fit_mat_2f


## CFI plot

cfi_dat_1f = fit_mat_1f[1, 1:n.cnts]
cfi_ord = order(as.numeric(cfi_dat_1f))
cfi_dat_1f = cfi_dat_1f[, cfi_ord]
names_cfi_1f = names(cfi_dat_1f)

cfi_dat_2f = fit_mat_2f[1, 1:n.cnts]
cfi_dat_2f = cfi_dat_2f[, cfi_ord]
names_cfi_2f = names(cfi_dat_2f)

pdf("fig_1a.pdf", height=5, width=3)
par(mfrow=c(1,1), mar=c(2.5, 5, 0.5, 1), tcl=-0.2)
plot(x=as.numeric(cfi_dat_1f), y=1:n.cnts, xlim=c(0.7, 1), type="n", axes=FALSE, 
     xlab="", ylab="")
axis(side=1, at=seq(0.7, 1, by=0.1), mgp=c(1, 0.2, 0), cex.axis=0.8)
axis(side=2, at=1:n.cnts, labels=names_cfi_1f, las=1, mgp=c(1, 0.4, 0), cex.axis=0.8)
mtext(side=1, line=1.2, text="CFI")
abline(h=1:n.cnts, lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=seq(0.7, 1, by=0.1), lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=0.90, lty=1, col=rgb(0.5, 0.5, 0.5))
points(x=as.numeric(cfi_dat_1f), y=1:n.cnts+0.2, type="p", pch=19, col=rgb(0, 0.4, 0.4, 1), 
       cex=1.2)
points(x=as.numeric(cfi_dat_2f), y=1:n.cnts-0.2, type="p", pch=19, col=rgb(1, 0.55, 0, 1), 
       cex=1.2)
segments(x0=0, x1=as.numeric(cfi_dat_1f), y0=1:n.cnts+0.2, lwd=1.5, col=rgb(0, 0.4, 0.4, 1))
segments(x0=0, x1=as.numeric(cfi_dat_2f), y0=1:n.cnts-0.2, lwd=1.5, col=rgb(1, 0.55, 0, 1))
box()
dev.off()


## RMSEA plot

rmsea_dat_1f = fit_mat_1f[3, 1:n.cnts]
rmsea_ord = order(as.numeric(rmsea_dat_1f))
rmsea_dat_1f = rmsea_dat_1f[, rmsea_ord]
names_rmsea_1f = names(rmsea_dat_1f)

rmsea_dat_2f = fit_mat_2f[3, 1:n.cnts]
rmsea_dat_2f = rmsea_dat_2f[, rmsea_ord]
names_rmsea_2f = names(rmsea_dat_2f)

pdf("fig_1b.pdf", height=5, width=3)
par(mfrow=c(1,1), mar=c(2.5, 5, 0.5, 1), tcl=-0.2)
plot(x=as.numeric(rmsea_dat_1f), y=n.cnts:1, xlim=c(0, 0.15), type="n", axes=FALSE, 
     xlab="", ylab="")
axis(side=1, at=seq(0, 0.4, by=0.05), mgp=c(1, 0.2, 0), cex.axis=0.8)
axis(side=2, at=n.cnts:1, labels=names_rmsea_1f, las=1, mgp=c(1, 0.4, 0), cex.axis=0.8)
mtext(side=1, line=1.2, text="RMSEA")
abline(h=1:n.cnts, lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=seq(0, 0.15, by=0.05), lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=0.08, lty=1, col=rgb(0.5, 0.5, 0.5))
points(x=as.numeric(rmsea_dat_1f), y=n.cnts:1+0.2, type="p", pch=19, 
       col=rgb(0, 0.4, 0.4, 1), cex=1.2)
points(x=as.numeric(rmsea_dat_2f), y=n.cnts:1-0.2, type="p", pch=19, 
       col=rgb(1, 0.55, 0, 1), cex=1.2)
segments(x0=1, x1=as.numeric(rmsea_dat_1f), y0=n.cnts:1+0.2, lwd=1.5, 
         col=rgb(0, 0.4, 0.4, 1))
segments(x0=1, x1=as.numeric(rmsea_dat_2f), y0=n.cnts:1-0.2, lwd=1.5, 
         col=rgb(1, 0.55, 0, 1))
legend("topright", legend=c("1-factor", "2-factor"), pch=19, bg=rgb(1, 1, 1, 0.75), 
       col=c(rgb(0, 0.4, 0.4, 1), rgb(1, 0.55, 0, 1)), box.col="white", 
       x.intersp=0.75, cex=0.8, pt.cex=1.2)
box()
dev.off()


## SRMR plot

srmr_dat_1f = fit_mat_1f[5, 1:n.cnts]
srmr_ord = order(as.numeric(srmr_dat_1f))
srmr_dat_1f = srmr_dat_1f[, srmr_ord]
names_srmr_1f = names(srmr_dat_1f)

srmr_dat_2f = fit_mat_2f[5, 1:n.cnts]
srmr_dat_2f = srmr_dat_2f[, srmr_ord]
names_srmr_2f = names(srmr_dat_2f)

pdf("fig_1c.pdf", height=5, width=3)
par(mfrow=c(1,1), mar=c(2.5, 5, 0.5, 1), tcl=-0.2)
plot(x=as.numeric(srmr_dat_1f), y=n.cnts:1, xlim=c(0, 0.15), type="n", 
     axes=FALSE, xlab="", ylab="")
axis(side=1, at=seq(0, 0.25, by=0.05), mgp=c(1, 0.2, 0), cex.axis=0.8)
axis(side=2, at=n.cnts:1, labels=names_srmr_1f, las=1, mgp=c(1, 0.4, 0), cex.axis=0.8)
mtext(side=1, line=1.2, text="SRMR")
abline(h=1:n.cnts, lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=seq(0, 0.15, by=0.05), lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=0.08, lty=1, col=rgb(0.5, 0.5, 0.5))
points(x=as.numeric(srmr_dat_1f), y=n.cnts:1+0.2, type="p", pch=19, 
       col=rgb(0, 0.4, 0.4, 1), cex=1.2)
points(x=as.numeric(srmr_dat_2f), y=n.cnts:1-0.2, type="p", pch=19, 
       col=rgb(1, 0.55, 0, 1), cex=1.2)
segments(x0=1, x1=as.numeric(srmr_dat_1f), y0=n.cnts:1+0.2, lwd=1.5, 
         col=rgb(0, 0.4, 0.4, 1))
segments(x0=1, x1=as.numeric(srmr_dat_2f), y0=n.cnts:1-0.2, lwd=1.5, 
         col=rgb(1, 0.55, 0, 1))
box()
dev.off()




#### Plot correlations between 17-item fit metrics and auxiliary variables

# read country covariate data
cov_dat = read.csv("country covariates.csv")

# rearrange datasets
fit_mat_t = as.data.frame(t(fit_mat_1f))
fit_df = merge(cov_dat, fit_mat_t, by.x = "country", by.y = 0)

# set up x and y vars for plot
x_vars = c("log(fit_df$gdp_cp_ppp)", "fit_df$ldi", "fit_df$eth_frac", "fit_df$internet_pc")
y_vars = c("fit_df$cfi", "fit_df$srmr", "fit_df$rmsea")
x_labs = c("log GDP per capita", "Liberal democracy index", "Ethnic fractionalization index", "Internet usage (%)")
y_labs = c("CFI", "SRMR", "RMSEA")

## Country covariate vs model fit plots

pdf("fig_s2.pdf", height = 8, width = 6)
par(mfrow = c(4, 3), mar = c(2.2, 2.0, 0.2, 0.5), tcl = -0.2, cex = 0.9, 
    cex.axis = 0.65)

for(i in 1:length(x_vars)) {
  for(j in 1:length(y_vars)) {
    x = as.numeric(eval(parse(text = x_vars[i])))
    y = as.numeric(eval(parse(text = y_vars[j])))
    plot(x, y, pch = 16, col = rgb(0, 0, 0, 0.5), 
         xlab = "", ylab = "", axes = FALSE)
    lines(lowess(x, y), lwd = 2, col = rgb(1, 0.55, 0, 0.9))
    abline(lm(y ~ x), lwd = 2, lty = 2, col = rgb(0, 0.4, 0.4, 0.9))
    axis(side = 1, labels = TRUE, mgp = c(1, -0.1, 0))
    axis(side = 2, labels = TRUE, mgp = c(1, 0.2, 0))
    mtext(side = 1, line = 0.7, text = x_labs[i], cex = 0.7)
    mtext(side = 2, line = 0.9, text = y_labs[j], cex = 0.7)
    box()
  }
}

dev.off()




#### 7-item CFA fit plot

fit_list_7 = list()

# Store data from each country's results in a list
for (cnt in cnts) {
  file_path = paste0(cnt, "/cfa_7it_ord_std_fit.csv")
  fit_dat_7 = read.csv(file_path, row.names = 1)
  fit_list_7[[cnt]] = fit_dat_7
  
}

# Combine the data into a matrix
fit_mat_7 = do.call(cbind, fit_list_7)
names(fit_mat_7) = cnts
fit_mat_7


## CFI plot

cfi_dat_7 = fit_mat_7[1, 1:n.cnts]
cfi_dat_7 = cfi_dat_7[, order(as.numeric(cfi_dat_7))]
names_cfi = names(cfi_dat_7)

pdf("fig_3a.pdf", height=5, width=3)
par(mfrow=c(1,1), mar=c(2.5, 5, 0.5, 1), tcl=-0.2)
plot(x=as.numeric(cfi_dat_7), y=1:n.cnts, xlim=c(0.85, 1), type="n", axes=FALSE, xlab="", ylab="")
axis(side=1, at=seq(0.75, 1, by=0.05), mgp=c(1, 0.2, 0), cex.axis=0.8)
axis(side=2, at=1:n.cnts, labels=names_cfi, las=1, mgp=c(1, 0.4, 0), cex.axis=0.8)
mtext(side=1, line=1.2, text="CFI")
abline(h=1:n.cnts, lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=seq(0.85, 1, by=0.05), lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=0.90, lty=1, col=rgb(0.5, 0.5, 0.5))
points(x=as.numeric(cfi_dat_7), y=1:n.cnts, type="p", pch=19, col=rgb(0, 0.4, 0.4, 1), cex=1.2)
segments(x0=0, x1=as.numeric(cfi_dat_7), y0=1:n.cnts, lwd=1.5, col=rgb(0, 0.4, 0.4, 1))
box()
dev.off()


## RMSEA plot

rmsea_dat_7 = fit_mat_7[3, 1:n.cnts]
rmsea_dat_7 = rmsea_dat_7[, order(as.numeric(rmsea_dat_7))]
names_rmsea = names(rmsea_dat_7)

pdf("fig_3b.pdf", height=5, width=3)
par(mfrow=c(1,1), mar=c(2.5, 5, 0.5, 1), tcl=-0.2)
plot(x=as.numeric(rmsea_dat_7), y=n.cnts:1, xlim=c(0, 0.11), type="n", axes=FALSE, xlab="", ylab="")
axis(side=1, at=seq(0, 0.11, by=0.02), mgp=c(1, 0.2, 0), cex.axis=0.8)
axis(side=2, at=n.cnts:1, labels=names_rmsea, las=1, mgp=c(1, 0.4, 0), cex.axis=0.8)
mtext(side=1, line=1.2, text="RMSEA")
abline(h=1:n.cnts, lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=seq(0, 0.11, by=0.02), lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=0.08, lty=1, col=rgb(0.5, 0.5, 0.5))
points(x=as.numeric(rmsea_dat_7), y=n.cnts:1, type="p", pch=19, col=rgb(0, 0.4, 0.4, 1), cex=1.2)
segments(x0=1, x1=as.numeric(rmsea_dat_7), y0=n.cnts:1, lwd=1.5, col=rgb(0, 0.4, 0.4, 1))
box()
dev.off()


## SRMR plot

srmr_dat_7 = fit_mat_7[5, 1:n.cnts]
srmr_dat_7 = srmr_dat_7[, order(as.numeric(srmr_dat_7))]
names_srmr = names(srmr_dat_7)

pdf("fig_3c.pdf", height=5, width=3)
par(mfrow=c(1,1), mar=c(2.5, 5, 0.5, 1), tcl=-0.2)
plot(x=as.numeric(srmr_dat_7), y=n.cnts:1, xlim=c(0, 0.10), type="n", axes=FALSE, xlab="", ylab="")
axis(side=1, at=seq(0, 0.10, by=0.02), mgp=c(1, 0.2, 0), cex.axis=0.8)
axis(side=2, at=n.cnts:1, labels=names_srmr, las=1, mgp=c(1, 0.4, 0), cex.axis=0.8)
mtext(side=1, line=1.2, text="SRMR")
abline(h=1:n.cnts, lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=seq(0, 0.10, by=0.02), lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=0.08, lty=1, col=rgb(0.5, 0.5, 0.5))
points(x=as.numeric(srmr_dat_7), y=n.cnts:1, type="p", pch=19, col=rgb(0, 0.4, 0.4, 1), cex=1.2)
segments(x0=1, x1=as.numeric(srmr_dat_7), y0=n.cnts:1, lwd=1.5, col=rgb(0, 0.4, 0.4, 1))
box()
dev.off()




#### 7-item cronbach's alpha

alph_vec_7 = vector("numeric", length = n.cnts)

# Store data from each country's results in a list
for (i in 1:length(cnts)) {
  file_path = paste0(cnts[i], "/sd7_ord_alpha.csv")
  alph_dat_7 = read.csv(file_path, header=TRUE)
  alph_vec_7[i] = alph_dat_7[2][[1]]
  
}

# add country names and reorder
names(alph_vec_7) = cnts
alph_dat_7 = alph_vec_7[order(as.numeric(alph_vec_7))]
names_alph = names(alph_dat_7)


## Cronbach's alpha plot

pdf("fig_4b.pdf", height=5, width=3)
par(mfrow=c(1,1), mar=c(2.5, 5, 0.5, 1), tcl=-0.2)
plot(x=as.numeric(alph_dat_7), y=1:n.cnts, xlim=c(0.5, 1), type="n", axes=FALSE, xlab="", ylab="")
axis(side=1, at=seq(0.5, 1, by=0.1), mgp=c(1, 0.2, 0), cex.axis=0.8)
axis(side=2, at=1:n.cnts, labels=names_alph, las=1, mgp=c(1, 0.4, 0), cex.axis=0.8)
mtext(side=1, line=1.2, text="Cronbach's alpha")
abline(h=1:n.cnts, lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=seq(0.5, 1, by=0.1), lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=0.70, lty=1, col=rgb(0.5, 0.5, 0.5))
points(x=as.numeric(alph_dat_7), y=1:n.cnts, type="p", pch=19, col=rgb(0, 0.4, 0.4, 1), cex=1.2)
segments(x0=0, x1=as.numeric(alph_dat_7), y0=1:n.cnts, lwd=1.5, col=rgb(0, 0.4, 0.4, 1))
box()
dev.off()



#### Correlations with criterion variables

# Initialize lists to store the data
cormat_list = list()
cormat_list_linz = list()
cormat_list_church = list()
cormat_list_imp = list()
cormat_list_strong = list()

# Store data from each country's results in a list
for (cnt in cnts) {
  file_path = paste0(cnt, "/sd7_crit_cormat.csv")
  
  # Check if the file exists
  if (file.exists(file_path)) {
    cormat = read.csv(file_path, row.names = 1)
    
    # Extract correlations for 7-item scale
    cormat_list[[cnt]] = cormat[, "SUPDEM_7IT", drop = FALSE]
    
    # Extract correlations for other sup dem items
    if ("LINZ_SUPDEM" %in% colnames(cormat)) {
      cormat_list_linz[[cnt]] = cormat[, "LINZ_SUPDEM", drop = FALSE]
    }
    if ("CHURCH_SUPDEM" %in% colnames(cormat)) {
      cormat_list_church[[cnt]] = cormat[, "CHURCH_SUPDEM", drop = FALSE]
    }
    if ("IMPORT_DEMOC" %in% colnames(cormat)) {
      cormat_list_imp[[cnt]] = cormat[, "IMPORT_DEMOC", drop = FALSE]
    }
    if ("STRONG_LEAD" %in% colnames(cormat)) {
      cormat_list_strong[[cnt]] = cormat[, "STRONG_LEAD", drop = FALSE]
    }
  }
}
    
# Set up an empty matrix with the specified row names
row_names = c("SUPDEM_7IT", "LINZ_SUPDEM", "CHURCH_SUPDEM", "STRONG_LEAD", "IMPORT_DEMOC",
              "SATIS_DEM", "GOVT_APPROV", "PRESIDENT_APPROV", "EXEC_APPROV", "LR_IDEOL",
              "LC_IDEOL", "POP_ATT", "POL_TRUST")
cormat_comb = matrix(nrow = length(row_names), ncol = length(cnts))
rownames(cormat_comb) = row_names
colnames(cormat_comb) = cnts

# Loop through the countries and add their data to the matrix
for (i in seq_along(cnts)) {
  cnt = cnts[i]
  
  # Check if data exists for this country
  if (!is.null(cormat_list[[cnt]])) {
    # Add the data to the matrix
    cormat_comb[, i] = cormat_list[[cnt]][row_names, , drop = TRUE]
  }
}

# Recode some variables
cormat_comb["EXEC_APPROV",] = ifelse(!(is.na(cormat_comb["GOVT_APPROV",])), 
                                     cormat_comb["GOVT_APPROV",], cormat_comb["PRESIDENT_APPROV",])
cormat_comb["LR_IDEOL",] = ifelse(is.na(cormat_comb["LR_IDEOL",]), 
                                  cormat_comb["LC_IDEOL",], cormat_comb["LR_IDEOL",])

# Reorder, drop missing countries, and reshape
cormat_comb = cormat_comb[c("LINZ_SUPDEM", "CHURCH_SUPDEM", "IMPORT_DEMOC", "STRONG_LEAD", "SATIS_DEM", 
                             "POL_TRUST", "POP_ATT", "EXEC_APPROV", "LR_IDEOL"), ]
cormat_comb = cormat_comb[, colSums(!is.na(cormat_comb)) > 0]
cormat_comb = cormat_comb[nrow(cormat_comb):1, ]
cormat_long = reshape2::melt(cormat_comb)

# Set up plot optioons
plot_col = brewer.pal(5, "Set2")
plot_col = c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00")
plot_pch = c(15, 16, 18)
plot_comb = expand.grid(plot_col, plot_pch, stringsAsFactors = FALSE)
plot_comb$Var3 = ifelse(plot_comb$Var2 == 15, 0.9, 
                         ifelse(plot_comb$Var2 == 16, 1, 1.1))

# axis labels
y_labs = c("Conservative\nideology", "Executive\napproval", "Populist\nattitudes", "Political\ntrust",
           "Satisfaction\nwith democracy", "Supports\nstrong leader", "Democracy is\nimportant", 
           "Democracy\nis best", "Democracy\npreferred")

# save correlation matrix
cormat_out = round(t(cormat_comb), 2)
colnames(cormat_out) = y_labs
write.csv(cormat_out, file = "table_s1.csv", row.names = TRUE)


## Country name dotplot

cnts_code = countrycode(colnames(cormat_comb), origin="country.name", destination="iso2c")

pdf("fig_4a.pdf", width=6, height=5)
par(mfrow = c(1, 1), mar = c(2.5, 5.5, 0.5, 0.5), cex = 1, las = 0, tcl = -0.2)
plot(NA, xlim = c(-0.8, 0.8), ylim = c(0.5, nrow(cormat_comb) + 0.5), xaxt = "n", yaxt = "n", bty = "n",
     xlab = "", ylab = "")
axis(1, at = seq(-1, +1, by = .2), mgp = c(1, 0.2, 0), cex.axis = 0.8)
axis(2, at = seq(1, nrow(cormat_comb)), mgp = c(1, 0.4, 0), labels = y_labs, las = 2, 
     cex.axis = 0.8, tcl = 0)
abline(v = seq(-1, +1, by = .2), lty = 3, col = rgb(0.5, 0.5, 0.5, 0.75))
abline(v = 0, lty = 1, col = rgb(0.5, 0.5, 0.5, 1))
abline(h = seq(0.5, nrow(cormat_comb) + 0.5), lty = 3, col = rgb(0.5, 0.5, 0.5, 0.75))
#abline(h = 1:nrow(cormat_comb), lty = 3, col = rgb(0.5, 0.5, 0.5, 0.75))
mtext(side = 1, line = 1.2, "Correlations with concise support for democracy scale", cex = 1)
for (i in seq_along(colnames(cormat_comb))) {
  country_data = subset(cormat_long, Var2 == colnames(cormat_comb)[i])
  text(x = country_data$value, y = jitter(as.numeric(country_data$Var1), factor = 1.5),  
       labels = cnts_code[i], col = rgb(0, 0, 0, 1), cex=0.65)
  # text(x = country_data$value, y = as.numeric(country_data$Var1) + (i - (ncol(cormat_comb) + 1) / 2) * 0.05, 
  #      labels = cnts_code[i], col = rgb(0, 0, 0, 1), cex=0.7)
}
box()
dev.off()



#### Dotplots for all measures of support for democracy

## wrangle data for the linz item

cormat_comb_linz = matrix(nrow = length(row_names), ncol = length(cnts))
rownames(cormat_comb_linz) = row_names
colnames(cormat_comb_linz) = cnts

for (i in seq_along(cnts)) {
  cnt = cnts[i]
  if (!is.null(cormat_list_linz[[cnt]])) {
    cormat_comb_linz[, i] = cormat_list_linz[[cnt]][row_names, , drop = TRUE]
  }
}

cormat_comb_linz["EXEC_APPROV",] = ifelse(!(is.na(cormat_comb_linz["GOVT_APPROV",])), 
                                          cormat_comb_linz["GOVT_APPROV",], 
                                          cormat_comb_linz["PRESIDENT_APPROV",])
cormat_comb_linz["LR_IDEOL",] = ifelse(is.na(cormat_comb_linz["LR_IDEOL",]), 
                                       cormat_comb_linz["LC_IDEOL",], 
                                       cormat_comb_linz["LR_IDEOL",])
cormat_comb_linz = cormat_comb_linz[
  c("LINZ_SUPDEM", "CHURCH_SUPDEM", "IMPORT_DEMOC", "STRONG_LEAD", "SATIS_DEM", 
    "POL_TRUST", "POP_ATT", "EXEC_APPROV", "LR_IDEOL"), ]

cormat_comb_linz = cormat_comb_linz[, colSums(!is.na(cormat_comb_linz)) > 0]
cormat_comb_linz = cormat_comb_linz[nrow(cormat_comb_linz):1, ]
cormat_long_linz = reshape2::melt(cormat_comb_linz)


## wrangle data for the churchill item

cormat_comb_church = matrix(nrow = length(row_names), ncol = length(cnts))
rownames(cormat_comb_church) = row_names
colnames(cormat_comb_church) = cnts

for (i in seq_along(cnts)) {
  cnt = cnts[i]
  if (!is.null(cormat_list_church[[cnt]])) {
    cormat_comb_church[, i] = cormat_list_church[[cnt]][row_names, , drop = TRUE]
  }
}

cormat_comb_church["EXEC_APPROV",] = ifelse(!(is.na(cormat_comb_church["GOVT_APPROV",])), 
                                          cormat_comb_church["GOVT_APPROV",], 
                                          cormat_comb_church["PRESIDENT_APPROV",])
cormat_comb_church["LR_IDEOL",] = ifelse(is.na(cormat_comb_church["LR_IDEOL",]), 
                                       cormat_comb_church["LC_IDEOL",], 
                                       cormat_comb_church["LR_IDEOL",])
cormat_comb_church = cormat_comb_church[
  c("LINZ_SUPDEM", "CHURCH_SUPDEM", "IMPORT_DEMOC", "STRONG_LEAD", "SATIS_DEM", 
    "POL_TRUST", "POP_ATT", "EXEC_APPROV", "LR_IDEOL"), ]

cormat_comb_church = cormat_comb_church[, colSums(!is.na(cormat_comb_church)) > 0]
cormat_comb_church = cormat_comb_church[nrow(cormat_comb_church):1, ]
cormat_long_church = reshape2::melt(cormat_comb_church)


## wrangle data for the importance item

cormat_comb_imp = matrix(nrow = length(row_names), ncol = length(cnts))
rownames(cormat_comb_imp) = row_names
colnames(cormat_comb_imp) = cnts

for (i in seq_along(cnts)) {
  cnt = cnts[i]
  if (!is.null(cormat_list_imp[[cnt]])) {
    cormat_comb_imp[, i] = cormat_list_imp[[cnt]][row_names, , drop = TRUE]
  }
}

cormat_comb_imp["EXEC_APPROV",] = ifelse(!(is.na(cormat_comb_imp["GOVT_APPROV",])), 
                                          cormat_comb_imp["GOVT_APPROV",], 
                                          cormat_comb_imp["PRESIDENT_APPROV",])
cormat_comb_imp["LR_IDEOL",] = ifelse(is.na(cormat_comb_imp["LR_IDEOL",]), 
                                       cormat_comb_imp["LC_IDEOL",], 
                                       cormat_comb_imp["LR_IDEOL",])
cormat_comb_imp = cormat_comb_imp[
  c("LINZ_SUPDEM", "CHURCH_SUPDEM", "IMPORT_DEMOC", "STRONG_LEAD", "SATIS_DEM", 
    "POL_TRUST", "POP_ATT", "EXEC_APPROV", "LR_IDEOL"), ]

cormat_comb_imp = cormat_comb_imp[, colSums(!is.na(cormat_comb_imp)) > 0]
cormat_comb_imp = cormat_comb_imp[nrow(cormat_comb_imp):1, ]
cormat_long_imp = reshape2::melt(cormat_comb_imp)


## wrangle data for the strong leader item

cormat_comb_strong = matrix(nrow = length(row_names), ncol = length(cnts))
rownames(cormat_comb_strong) = row_names
colnames(cormat_comb_strong) = cnts

for (i in seq_along(cnts)) {
  cnt = cnts[i]
  if (!is.null(cormat_list_strong[[cnt]])) {
    cormat_comb_strong[, i] = cormat_list_strong[[cnt]][row_names, , drop = TRUE]
  }
}

cormat_comb_strong["EXEC_APPROV",] = ifelse(!(is.na(cormat_comb_strong["GOVT_APPROV",])), 
                                          cormat_comb_strong["GOVT_APPROV",], 
                                          cormat_comb_strong["PRESIDENT_APPROV",])
cormat_comb_strong["LR_IDEOL",] = ifelse(is.na(cormat_comb_strong["LR_IDEOL",]), 
                                       cormat_comb_strong["LC_IDEOL",], 
                                       cormat_comb_strong["LR_IDEOL",])
cormat_comb_strong = cormat_comb_strong[
  c("LINZ_SUPDEM", "CHURCH_SUPDEM", "IMPORT_DEMOC", "STRONG_LEAD", "SATIS_DEM", 
    "POL_TRUST", "POP_ATT", "EXEC_APPROV", "LR_IDEOL"), ]

cormat_comb_strong = cormat_comb_strong[, colSums(!is.na(cormat_comb_strong)) > 0]
cormat_comb_strong = cormat_comb_strong[nrow(cormat_comb_strong):1, ]
cormat_long_strong = reshape2::melt(cormat_comb_strong)


## Plot all 5

# axis labels
y_labs = c("Conservative ideol.", "Executive approval", "Populist attitudes", "Political trust",
           "Satisf. w/ democracy", "Supp. strong leader", "Democ. is important", 
           "Democracy is best", "Democracy preferred")


pdf("fig_s3.pdf", height = 8, width = 8)

par(mfrow = c(3, 2), mar = c(1.5, 6, 2.5, 0.5), cex = 1, las = 0, tcl = -0.2)

# 7-item plot
plot(NA, xlim = c(-0.8, 0.8), ylim = c(0.5, nrow(cormat_comb) + 0.5), xaxt = "n", yaxt = "n", bty = "n",
     xlab = "", ylab = "")
axis(1, at = seq(-1, +1, by = .2), mgp = c(1, 0.2, 0), cex.axis = 0.65)
axis(2, at = seq(1, nrow(cormat_comb)), mgp = c(1, 0.4, 0), labels = y_labs, las = 2, 
     cex.axis = 0.65, tcl = 0)
abline(v = seq(-1, +1, by = .2), lty = 3, col = rgb(0.5, 0.5, 0.5, 0.75))
abline(v = 0, lty = 1, col = rgb(0.5, 0.5, 0.5, 1))
abline(h = seq(0.5, nrow(cormat_comb) + 0.5), lty = 3, col = rgb(0.5, 0.5, 0.5, 0.75))
mtext(side = 3, line = 0.3, "Correlations with\nconcise scale", cex = 0.9)
for (i in seq_along(colnames(cormat_comb))) {
  country_data = subset(cormat_long, Var2 == colnames(cormat_comb)[i])
  points(country_data$value, as.numeric(country_data$Var1), 
         col = rgb(0, 0, 0, 0.6), pch = 16, cex = 1)
}
box()

# linz plot
plot(NA, xlim = c(-0.8, 0.8), ylim = c(0.5, nrow(cormat_comb_linz) + 0.5), xaxt = "n", 
     yaxt = "n", bty = "n", xlab = "", ylab = "")
axis(1, at = seq(-1, +1, by = .2), mgp = c(1, 0.2, 0), cex.axis = 0.65)
axis(2, at = seq(1, nrow(cormat_comb_linz)), mgp = c(1, 0.4, 0), labels = y_labs, las = 2, 
     cex.axis = 0.65, tcl = 0)
abline(v = seq(-1, +1, by = .2), lty = 3, col = rgb(0.5, 0.5, 0.5, 0.75))
abline(v = 0, lty = 1, col = rgb(0.5, 0.5, 0.5, 1))
abline(h = seq(0.5, nrow(cormat_comb_linz) + 0.5), lty = 3, col = rgb(0.5, 0.5, 0.5, 0.75))
mtext(side = 3, line = 0.3, "Correlations with\n'democracy preferred' item", cex = 0.9)
for (i in seq_along(colnames(cormat_comb_linz))) {
  country_data = subset(cormat_long_linz, Var2 == colnames(cormat_comb_linz)[i])
  points(country_data$value, as.numeric(country_data$Var1), 
         col = rgb(0, 0, 0, 0.6), pch = 16, cex = 1)
}
box()

# church plot
plot(NA, xlim = c(-0.8, 0.8), ylim = c(0.5, nrow(cormat_comb_church) + 0.5), xaxt = "n", 
     yaxt = "n", bty = "n", xlab = "", ylab = "")
axis(1, at = seq(-1, +1, by = .2), mgp = c(1, 0.2, 0), cex.axis = 0.65)
axis(2, at = seq(1, nrow(cormat_comb_church)), mgp = c(1, 0.4, 0), labels = y_labs, las = 2, 
     cex.axis = 0.65, tcl = 0)
abline(v = seq(-1, +1, by = .2), lty = 3, col = rgb(0.5, 0.5, 0.5, 0.75))
abline(v = 0, lty = 1, col = rgb(0.5, 0.5, 0.5, 1))
abline(h = seq(0.5, nrow(cormat_comb_church) + 0.5), lty = 3, col = rgb(0.5, 0.5, 0.5, 0.75))
mtext(side = 3, line = 0.3, "Correlations with\n'democracy is best' item", cex = 0.9)
for (i in seq_along(colnames(cormat_comb_church))) {
  country_data = subset(cormat_long_church, Var2 == colnames(cormat_comb_church)[i])
  points(country_data$value, as.numeric(country_data$Var1), 
         col = rgb(0, 0, 0, 0.6), pch = 16, cex = 1)
}
box()

# importance plot
plot(NA, xlim = c(-0.8, 0.8), ylim = c(0.5, nrow(cormat_comb_imp) + 0.5), xaxt = "n", 
     yaxt = "n", bty = "n", xlab = "", ylab = "")
axis(1, at = seq(-1, +1, by = .2), mgp = c(1, 0.2, 0), cex.axis = 0.65)
axis(2, at = seq(1, nrow(cormat_comb_imp)), mgp = c(1, 0.4, 0), labels = y_labs, las = 2, 
     cex.axis = 0.65, tcl = 0)
abline(v = seq(-1, +1, by = .2), lty = 3, col = rgb(0.5, 0.5, 0.5, 0.75))
abline(v = 0, lty = 1, col = rgb(0.5, 0.5, 0.5, 1))
abline(h = seq(0.5, nrow(cormat_comb_imp) + 0.5), lty = 3, col = rgb(0.5, 0.5, 0.5, 0.75))
mtext(side = 3, line = 0.3, "Correlations with\n'democracy is important' item", cex = 0.9)
for (i in seq_along(colnames(cormat_comb_imp))) {
  country_data = subset(cormat_long_imp, Var2 == colnames(cormat_comb_imp)[i])
  points(country_data$value, as.numeric(country_data$Var1), 
         col = rgb(0, 0, 0, 0.6), pch = 16, cex = 1)
}
box()

# strong plot
plot(NA, xlim = c(-0.8, 0.8), ylim = c(0.5, nrow(cormat_comb_strong) + 0.5), xaxt = "n", 
     yaxt = "n", bty = "n", xlab = "", ylab = "")
axis(1, at = seq(-1, +1, by = .2), mgp = c(1, 0.2, 0), cex.axis = 0.65)
axis(2, at = seq(1, nrow(cormat_comb_strong)), mgp = c(1, 0.4, 0), labels = y_labs, las = 2, 
     cex.axis = 0.65, tcl = 0)
abline(v = seq(-1, +1, by = .2), lty = 3, col = rgb(0.5, 0.5, 0.5, 0.75))
abline(v = 0, lty = 1, col = rgb(0.5, 0.5, 0.5, 1))
abline(h = seq(0.5, nrow(cormat_comb_strong) + 0.5), lty = 3, col = rgb(0.5, 0.5, 0.5, 0.75))
mtext(side = 3, line = 0.3, "Correlations with\n'strong leader' item", cex = 0.9)
for (i in seq_along(colnames(cormat_comb_strong))) {
  country_data = subset(cormat_long_strong, Var2 == colnames(cormat_comb_strong)[i])
  points(country_data$value, as.numeric(country_data$Var1),  
         col = rgb(0, 0, 0, 0.6), pch = 16, cex = 1)
}
box()

dev.off()



#### Plot eigenvalue results

merged_eigen = data.frame(matrix(NA, ncol=length(cnts), nrow=10))

# loop through each dataset
for(i in 1:length(cnts)) {
  file = paste0("./", cnts[i], "/sd_eigen.csv")
  data = read.csv(file)
  merged_eigen[,i] = data[1:10, 2]
  colnames(merged_eigen)[i] = paste0(cnts[i], "_eigen")
}

# title-case country names
cnt_plot = str_to_title(cnts)
cnt_plot[n.cnts] = "USA"

# plot colours
plot_cols = unname(alphabet(n=n.cnts))
plot_cols[4] = "#740AFF"
plot_cols[8] = "#FFA405"
plot_cols[10] = "#FF0010"
plot_rgb = hex2RGB(plot_cols)
plot_rgba = adjust_transparency(plot_rgb, alpha=0.6)


## Scree plots

pdf("fig_s1.pdf", height=5, width=8)
par(mfrow=c(1,1), mar=c(2.5, 2.5, 0.5, 1), tcl=-0.2)
plot(x=1:10, y=merged_eigen[,1], ylim=c(0, 5.5), yaxs="i", type="n", axes=FALSE, xlab="", ylab="")
axis(side=1, at=1:10, mgp=c(1, 0.2, 0), cex.axis=0.8)
axis(side=2, at=seq(0, 5, by=1), las=1, mgp=c(1, 0.4, 0), cex.axis=0.8)
mtext(side=1, line=1.2, text="Principal components")
mtext(side=2, line=1.4, text="Eigenvalues")
abline(v=1:10, lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(h=seq(0, 5, by=1), lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(h=1, lty=1, col=rgb(0.5, 0.5, 0.5))
for(i in 1:n.cnts) {
  points(x=1:10, y=merged_eigen[, i], type="b", pch=19, col=plot_rgba[i], lwd=1.5)
  points(x=1:10, y=merged_eigen[, i], type="p", pch=1, col=rgb(0.5, 0.5, 0.5, 1))
}
legend(x=9, y=5.4, legend=cnt_plot, col=plot_rgba, pch=19, bg="white", box.col="white", cex=0.9,
       x.intersp=0.6, y.intersp=0.9)
box()
dev.off()




#### Plot CFA loadings for weighted, continuous CFAs

n.cnts = length(cnts)

# define an empty list to store the tables
CFA_tables = list()

# loop over the country names
for (i in 1:length(cnts)) {
  
  # construct the file path
  file = paste0("./", cnts[i], "/cfa_1fac_std_table_w.html")
  
  # read the table using readHTMLTable()
  tab_temp = readHTMLTable(file, which = 1)
  
  # extract the columns V1 and V2 and rows 3 to 20
  tab_temp = tab_temp[3:20, c("V1", "V2")]
  
  # rename the column V2 with the country name
  colnames(tab_temp)[2] = cnts[i]
  
  # add the table to the list
  CFA_tables[[i]] = tab_temp
  
}

# merge the tables using V1 as the common key
CFA_merg = Reduce(function(x, y) merge(x, y, by = "V1", all = TRUE), CFA_tables)

# rearrange rows and clean up
row_lab = c("FREXP1", "FREXP2", "FRASSC1", "FRASSC2", "FRASSC3", "UNISUFF1", "UNISUFF2", "DECELEC1", 
            "DECELEC2", "FRELECT1", "FRELECT2", "JUDCNSTR1", "JUDCNSTR2", "LEGCNSTR1", "LEGCNSTR2", 
            "EQLAW1", "EQLAW2")
CFA_merg$V1 = factor(CFA_merg$V1, levels = row_lab)
CFA_merg = arrange(CFA_merg, V1)
CFA_merg = filter(CFA_merg, V1 != "")

# print the result
print(CFA_merg)

# remove the row names column
rownames(CFA_merg) = CFA_merg[,1]
CFA_merg[,1] = NULL


## Plot loadings

# country code
cnts1 = cnts
cnts_code = countrycode(cnts1, origin="country.name", destination="iso2c")

# title-case country names
cnt_plot = str_to_title(cnts1)
cnt_plot[n.cnts] = "USA"

# more meaningful labels
it_labs = c("FreeExp1(r)", "FreeExp2", "FreeAssc1", "FreeAssc2(r)", "FreeAssc3", "UniSuff1", "UniSuff2(r)", "ElecDecMk1", 
            "ElecDecMk2(r)", "FFElect1(r)", "FFElect2", "JudCnstr1(r)", "JudCnstr2", "LegCnstr1", "LegCnstr2(r)", "EqLaw1", 
            "EqLaw2(r)")


## Plot with component labels

comp_labs = c("Freedom of expression", "Freedom of association", "Universal suffrage", 
              "Key decision-makers elected", "Free and fair elections", "Judicial constraints",
              "Legislative constraints", "Equality before the law")

cfa_dat_2 = matrix(NA, nrow=25, ncol=19)
lab_inds = c(1, 4, 8, 11, 14, 17, 20, 23)
dat_inds = c(2, 3, 5, 6, 7, 9, 10, 12, 13, 15, 16, 18, 19, 21, 22, 24, 25)
comp_pos = c(25, 22, 18, 15, 12, 9, 6, 3)
lab_pos = c(24, 23, 21, 20, 19, 17, 16, 14, 13, 11, 10, 8, 7, 5, 4, 2, 1)
for(i in 1:n.cnts) {
  cfa_dat_2[lab_inds, i] = NA
  cfa_dat_2[dat_inds, i] = CFA_merg[, i] 
}
it_labs_2 = c("Criticize govt.(r)", "Free media", "One party", "Protest(r)", "Ban orgs", "Voter competence", 
              "Tolerance(r)", "Technocratic rule", "Unelected authority(r)", "Respect results(r)", 
              "Bend rules", "Judicial review(r)", "Judicial compliance", "Legisl. compliance", 
              "Legis. review(r)", "Law enforced", "Access justice(r)")
comp_headers = lab_inds

pdf("fig_s5.pdf", height=7, width=8)
par(mfrow=c(1,1), mar=c(2.5, 7, 0.5, 1), tcl=-0.2)
plot(x=cfa_dat_2[, 1], y=25:1, type="n", pch=16, col=rgb(0, 0, 0, 0.5), 
     xlim=c(-0.7, 1), xaxs="i", axes=FALSE, xlab="", ylab="")
axis(side=1, at=round(seq(-0.6, 0.8, by=0.2), 1), mgp=c(1, 0.2, 0), cex.axis=0.8)
axis(side=2, at=lab_pos, labels=it_labs_2, las=1, mgp=c(1, 0.4, 0), cex.axis=0.8)
mtext(side=1, line=1.2, text="Liberal democracy factor loading")
abline(h=lab_pos, lty=3, lwd=0.75, col=rgb(0.3, 0.3, 0.3, 0.75))
abline(v=seq(-0.6, 0.8, by=0.2), lty=3, lwd=0.75, col=rgb(0.3, 0.3, 0.3, 0.75))
abline(v=0, lty=1, lwd=0.75, col=rgb(0.3, 0.3, 0.3, 0.75))
abline(h=comp_pos, lty=1, lwd=25, col=rgb(0.7, 0.7, 0.7, 1))
text(x=0.2, y=comp_pos, adj=0.5, cex=0.9, labels=comp_labs)
for(i in 1:n.cnts){
  text(x=cfa_dat_2[, i], y=jitter(25:1, factor=1), labels=cnts_code[i], col=rgb(0, 0, 0, 1), cex=0.65)
}
box()
dev.off()



### 17-item CFA fit plot

## load 1 factor fits

# Store data from each country's results in a list
fit_list_1f = list()
for (cnt in cnts) {
  file_path = paste0(cnt, "/cfa_1fac_std_fit_w.csv")
  fit_dat_1f = read.csv(file_path, row.names = 1)
  fit_list_1f[[cnt]] = fit_dat_1f
}

# Combine the data from all files into a matrix
fit_mat_1f = do.call(cbind, fit_list_1f)
names(fit_mat_1f) = cnts
fit_mat_1f


## load 2 factor fits

# Store data from each country's results in a list
fit_list_2f = list()
for (cnt in cnts) {
  file_path = paste0(cnt, "/cfa_2fac_std_fit_w.csv")
  fit_dat_2f = read.csv(file_path, row.names = 1)
  fit_list_2f[[cnt]] = fit_dat_2f
}

# Combine the data from all files into a matrix
fit_mat_2f = do.call(cbind, fit_list_2f)
names(fit_mat_2f) = cnts
fit_mat_2f


## CFI plot

cfi_dat_1f = fit_mat_1f[1, 1:n.cnts]
cfi_ord = order(as.numeric(cfi_dat_1f))
cfi_dat_1f = cfi_dat_1f[, cfi_ord]
names_cfi_1f = names(cfi_dat_1f)

cfi_dat_2f = fit_mat_2f[1, 1:n.cnts]
cfi_dat_2f = cfi_dat_2f[, cfi_ord]
names_cfi_2f = names(cfi_dat_2f)

pdf("fig_s4a.pdf", height=5, width=3)
par(mfrow=c(1,1), mar=c(2.5, 5, 0.5, 1), tcl=-0.2)
plot(x=as.numeric(cfi_dat_1f), y=1:n.cnts, xlim=c(0.6, 1), type="n", axes=FALSE, 
     xlab="", ylab="")
axis(side=1, at=seq(0.6, 1, by=0.1), mgp=c(1, 0.2, 0), cex.axis=0.8)
axis(side=2, at=1:n.cnts, labels=names_cfi_1f, las=1, mgp=c(1, 0.4, 0), cex.axis=0.8)
mtext(side=1, line=1.2, text="CFI")
abline(h=1:n.cnts, lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=seq(0.5, 1, by=0.1), lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=0.90, lty=1, col=rgb(0.5, 0.5, 0.5))
points(x=as.numeric(cfi_dat_1f), y=1:n.cnts+0.2, type="p", pch=19, col=rgb(0, 0.4, 0.4, 1), 
       cex=1.2)
points(x=as.numeric(cfi_dat_2f), y=1:n.cnts-0.2, type="p", pch=19, col=rgb(1, 0.55, 0, 1), 
       cex=1.2)
segments(x0=0, x1=as.numeric(cfi_dat_1f), y0=1:n.cnts+0.2, lwd=1.5, col=rgb(0, 0.4, 0.4, 1))
segments(x0=0, x1=as.numeric(cfi_dat_2f), y0=1:n.cnts-0.2, lwd=1.5, col=rgb(1, 0.55, 0, 1))
box()
dev.off()

## RMSEA plot

rmsea_dat_1f = fit_mat_1f[3, 1:n.cnts]
rmsea_ord = order(as.numeric(rmsea_dat_1f))
rmsea_dat_1f = rmsea_dat_1f[, rmsea_ord]
names_rmsea_1f = names(rmsea_dat_1f)

rmsea_dat_2f = fit_mat_2f[3, 1:n.cnts]
rmsea_dat_2f = rmsea_dat_2f[, rmsea_ord]
names_rmsea_2f = names(rmsea_dat_2f)

pdf("fig_s4b.pdf", height=5, width=3)
par(mfrow=c(1,1), mar=c(2.5, 5, 0.5, 1), tcl=-0.2)
plot(x=as.numeric(rmsea_dat_1f), y=n.cnts:1, xlim=c(0, 0.15), type="n", axes=FALSE, 
     xlab="", ylab="")
axis(side=1, at=seq(0, 0.4, by=0.05), mgp=c(1, 0.2, 0), cex.axis=0.8)
axis(side=2, at=n.cnts:1, labels=names_rmsea_1f, las=1, mgp=c(1, 0.4, 0), cex.axis=0.8)
mtext(side=1, line=1.2, text="RMSEA")
abline(h=1:n.cnts, lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=seq(0, 0.15, by=0.05), lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=0.08, lty=1, col=rgb(0.5, 0.5, 0.5))
points(x=as.numeric(rmsea_dat_1f), y=n.cnts:1+0.2, type="p", pch=19, 
       col=rgb(0, 0.4, 0.4, 1), cex=1.2)
points(x=as.numeric(rmsea_dat_2f), y=n.cnts:1-0.2, type="p", pch=19, 
       col=rgb(1, 0.55, 0, 1), cex=1.2)
segments(x0=1, x1=as.numeric(rmsea_dat_1f), y0=n.cnts:1+0.2, lwd=1.5, 
         col=rgb(0, 0.4, 0.4, 1))
segments(x0=1, x1=as.numeric(rmsea_dat_2f), y0=n.cnts:1-0.2, lwd=1.5, 
         col=rgb(1, 0.55, 0, 1))
legend("topright", legend=c("1-factor", "2-factor"), pch=19, bg=rgb(1, 1, 1, 0.75), 
       col=c(rgb(0, 0.4, 0.4, 1), rgb(1, 0.55, 0, 1)), box.col="white", 
       x.intersp=0.75, cex=0.8, pt.cex=1.2)
box()
dev.off()

## SRMR plot

srmr_dat_1f = fit_mat_1f[5, 1:n.cnts]
srmr_ord = order(as.numeric(srmr_dat_1f))
srmr_dat_1f = srmr_dat_1f[, srmr_ord]
names_srmr_1f = names(srmr_dat_1f)

srmr_dat_2f = fit_mat_2f[5, 1:n.cnts]
srmr_dat_2f = srmr_dat_2f[, srmr_ord]
names_srmr_2f = names(srmr_dat_2f)

pdf("fig_s4c.pdf", height=5, width=3)
par(mfrow=c(1,1), mar=c(2.5, 5, 0.5, 1), tcl=-0.2)
plot(x=as.numeric(srmr_dat_1f), y=n.cnts:1, xlim=c(0, 0.15), type="n", 
     axes=FALSE, xlab="", ylab="")
axis(side=1, at=seq(0, 0.15, by=0.05), mgp=c(1, 0.2, 0), cex.axis=0.8)
axis(side=2, at=n.cnts:1, labels=names_srmr_1f, las=1, mgp=c(1, 0.4, 0), cex.axis=0.8)
mtext(side=1, line=1.2, text="SRMR")
abline(h=1:n.cnts, lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=seq(0, 0.25, by=0.05), lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=0.08, lty=1, col=rgb(0.5, 0.5, 0.5))
points(x=as.numeric(srmr_dat_1f), y=n.cnts:1+0.2, type="p", pch=19, 
       col=rgb(0, 0.4, 0.4, 1), cex=1.2)
points(x=as.numeric(srmr_dat_2f), y=n.cnts:1-0.2, type="p", pch=19, 
       col=rgb(1, 0.55, 0, 1), cex=1.2)
segments(x0=1, x1=as.numeric(srmr_dat_1f), y0=n.cnts:1+0.2, lwd=1.5, 
         col=rgb(0, 0.4, 0.4, 1))
segments(x0=1, x1=as.numeric(srmr_dat_2f), y0=n.cnts:1-0.2, lwd=1.5, 
         col=rgb(1, 0.55, 0, 1))
box()
dev.off()



#### 7-item CFA fit, weighted continuous models

n.cnts = length(cnts)
fit_list_7 = list()

# Store data from each country's results in a list
for (cnt in cnts) {
  file_path = paste0(cnt, "/cfa_7it_std_fit_w.csv")
  fit_dat_7 = read.csv(file_path, row.names = 1)
  fit_list_7[[cnt]] = fit_dat_7
  
}

# Combine the data into a matrix
fit_mat_7 = do.call(cbind, fit_list_7)
names(fit_mat_7) = cnts
fit_mat_7


## CFI plot

cfi_dat_7 = fit_mat_7[1, 1:n.cnts]
cfi_dat_7 = cfi_dat_7[, order(as.numeric(cfi_dat_7))]
names_cfi = names(cfi_dat_7)

pdf("fig_s6a.pdf", height=5, width=3)
par(mfrow=c(1,1), mar=c(2.5, 5, 0.5, 1), tcl=-0.2)
plot(x=as.numeric(cfi_dat_7), y=1:n.cnts, xlim=c(0.80, 1), type="n", axes=FALSE, xlab="", ylab="")
axis(side=1, at=seq(0.75, 1, by=0.05), mgp=c(1, 0.2, 0), cex.axis=0.8)
axis(side=2, at=1:n.cnts, labels=names_cfi, las=1, mgp=c(1, 0.4, 0), cex.axis=0.8)
mtext(side=1, line=1.2, text="CFI")
abline(h=1:n.cnts, lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=seq(0.80, 1, by=0.05), lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=0.90, lty=1, col=rgb(0.5, 0.5, 0.5))
points(x=as.numeric(cfi_dat_7), y=1:n.cnts, type="p", pch=19, col=rgb(0, 0.4, 0.4, 1), cex=1.2)
segments(x0=0, x1=as.numeric(cfi_dat_7), y0=1:n.cnts, lwd=1.5, col=rgb(0, 0.4, 0.4, 1))
box()
dev.off()


## RMSEA plot

rmsea_dat_7 = fit_mat_7[3, 1:n.cnts]
rmsea_dat_7 = rmsea_dat_7[, order(as.numeric(rmsea_dat_7))]
names_rmsea = names(rmsea_dat_7)

pdf("fig_s6b.pdf", height=5, width=3)
par(mfrow=c(1,1), mar=c(2.5, 5, 0.5, 1), tcl=-0.2)
plot(x=as.numeric(rmsea_dat_7), y=n.cnts:1, xlim=c(0, 0.11), type="n", axes=FALSE, xlab="", ylab="")
axis(side=1, at=seq(0, 0.11, by=0.02), mgp=c(1, 0.2, 0), cex.axis=0.8)
axis(side=2, at=n.cnts:1, labels=names_rmsea, las=1, mgp=c(1, 0.4, 0), cex.axis=0.8)
mtext(side=1, line=1.2, text="RMSEA")
abline(h=1:n.cnts, lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=seq(0, 0.11, by=0.02), lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=0.08, lty=1, col=rgb(0.5, 0.5, 0.5))
points(x=as.numeric(rmsea_dat_7), y=n.cnts:1, type="p", pch=19, col=rgb(0, 0.4, 0.4, 1), cex=1.2)
segments(x0=1, x1=as.numeric(rmsea_dat_7), y0=n.cnts:1, lwd=1.5, col=rgb(0, 0.4, 0.4, 1))
box()
dev.off()


## SRMR plot

srmr_dat_7 = fit_mat_7[5, 1:n.cnts]
srmr_dat_7 = srmr_dat_7[, order(as.numeric(srmr_dat_7))]
names_srmr = names(srmr_dat_7)

pdf("fig_s6c.pdf", height=5, width=3)
par(mfrow=c(1,1), mar=c(2.5, 5, 0.5, 1), tcl=-0.2)
plot(x=as.numeric(srmr_dat_7), y=n.cnts:1, xlim=c(0, 0.10), type="n", axes=FALSE, xlab="", ylab="")
axis(side=1, at=seq(0, 0.10, by=0.02), mgp=c(1, 0.2, 0), cex.axis=0.8)
axis(side=2, at=n.cnts:1, labels=names_srmr, las=1, mgp=c(1, 0.4, 0), cex.axis=0.8)
mtext(side=1, line=1.2, text="SRMR")
abline(h=1:n.cnts, lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=seq(0, 0.10, by=0.02), lty=3, col=rgb(0.5, 0.5, 0.5, 0.75))
abline(v=0.08, lty=1, col=rgb(0.5, 0.5, 0.5))
points(x=as.numeric(srmr_dat_7), y=n.cnts:1, type="p", pch=19, col=rgb(0, 0.4, 0.4, 1), cex=1.2)
segments(x0=1, x1=as.numeric(srmr_dat_7), y0=n.cnts:1, lwd=1.5, col=rgb(0, 0.4, 0.4, 1))
box()
dev.off()
